OISSTv2 Sea Surface Temperature Report
Adam A. Kemberling
Updated on: 2021-04-19
Apershing Gulf Of Maine Sea Surface Temperature Outlook
Region Extent
The spatial extent for Apershing Gulf Of Maine is displayed below. This bounding box is the same bounding box coordinates used to clip the OISST data when constructing the time series data from the array.
# testing regions
# params$region <- "apershing_gulf_of_maine"
# File paths for various extents based on params$region
region_paths <- get_timeseries_paths(region_group = "gmri_sst_focal_areas")
# Load the bounding box for Andy's GOM to show they align
poly_path <- region_paths[[params$region]][["shape_path"]]
region_extent <- st_read(poly_path, quiet = TRUE)
# Pull extents for the region for crop
crop_x <- st_bbox(region_extent)[c(1,3)]
crop_y <- st_bbox(region_extent)[c(2,4)]
# Zoom out for cpr extent
if(tolower(params$region) == "cpr gulf of maine"){
crop_x <- c(-70.875, -65.375)
crop_y <- c(40.375, 45.125)}
# one off shapes
# region_extent <- read_sf(paste0(res_path, ""))
# Full plot
ggplot() +
geom_sf(data = new_england, fill = "gray90") +
geom_sf(data = canada, fill = "gray90") +
geom_sf(data = greenland, fill = "gray90") +
geom_sf(data = region_extent,
color = gmri_cols("gmri blue"),
fill = gmri_cols("gmri blue"), alpha = 0.2, linetype = 2, size = 0.5) +
map_theme +
coord_sf(xlim = crop_x,
ylim = crop_y, expand = T) 
OISSTv2 Resources on Box
For many of the area we look at frequently here at GMRI, regional climatologies and anomalies have been processed for quick access to aid in other research endeavors.
These resources were processed using a handful of jupyter notebook work flows to take advantage of {xarray}.
The root folder for the following OISST products is ~/Box/RES_Data/OISST/oisst_mainstays. Resources can be accessed manually through the box folder or more directly using the gmRi R package.
Regional Timeseries
Area-specific time series are the most basic building block for relaying temporal trends. For any desired area (represented by a spatial polygon) we can generate a time series table of the mean sea surface temperature within that area for each day. Additionally, we can compare how observed temperatures correspond with the expected conditions based on a climatology using a specified reference period.
Accessing the Data
Pre-processed tables have been placed on box for quick access from within the RES_Data folder: RES_Data/OISST/oisst_mainstays/regional_timeseries.
A function to quickly access these time-lines has been added to the {gmRi} package as oisst_access_timeseries(). Additional functions exist to look up what timeseries are availablegmRi::get_region_names(). And an additional function was made to return both the path to the timeseries data and the path to the shapefile used to generate the timeseries using: gmRi::get_timeseries_paths().
# Use {gmRi} instead to load timeseries to tye up loose ends
timeseries_path <- region_paths[[params$region]][["timeseries_path"]]
region_timeseries <- read_csv(timeseries_path, col_types = cols())
# format dates
region_timeseries <- region_timeseries %>%
mutate(time = as.Date(time))
# Display Table of first 6 entries
head(region_timeseries) %>%
mutate(across(where(is.numeric), round, 2)) %>%
kable() %>%
kable_styling()| time | sst | modified_ordinal_day | sst_clim | clim_sd | sst_anom |
|---|---|---|---|---|---|
| 1981-09-01 | 15.78 | 245 | 16.43 | 2.27 | -0.65 |
| 1981-09-02 | 15.79 | 246 | 16.36 | 2.25 | -0.57 |
| 1981-09-03 | 15.49 | 247 | 16.30 | 2.23 | -0.81 |
| 1981-09-04 | 14.99 | 248 | 16.26 | 2.20 | -1.27 |
| 1981-09-05 | 14.84 | 249 | 16.15 | 2.18 | -1.30 |
| 1981-09-06 | 15.08 | 250 | 16.11 | 2.14 | -1.03 |
# march 1st sst
mar1 <- region_timeseries %>%
filter(modified_ordinal_day == 61) %>%
distinct(sst_clim) %>%
pull(sst_clim)Each of our Climatologies are currently set up to calculate daily averages on a modified year day, such that every March 1st and all days after fall on the same day, regardless of whether it is a leap year or not.
This preserves comparisons across calendar dates such-as: “The average temperature on march 1st is 4.2378197` for the reference period 1982 to 2011”
In these tables sst is the mean temperature observed for that date averaged across all cells within the area. sst_clim & clim_sd are the climate means and standard deviations for a 1982-2011 climatology. sst_anom is the daily observed minus the climate mean.
Warming Rates
Regional warming trends below were calculated using all the available data for complete years beginning with 1982 through the end of 2020.
Annual
# Summarize by year to return mean annual anomalies and variance
annual_summary <- region_timeseries %>%
mutate(year = year(time)) %>%
group_by(year) %>%
summarise(sst = mean(sst, na.rm = T),
sst_anom = mean(sst_anom, na.rm = T),
.groups = "keep") %>%
ungroup() %>%
mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")))
# # Global Temperature Anomaly Rates
global_anoms <- read_csv(paste0(oisst_path, "global_timeseries/global_anoms_1982to2011.csv"))
global_anoms <- mutate(global_anoms, year = year(time))
global_summary <- global_anoms %>%
group_by(year) %>%
summarise(sst_anom = mean(sst, na.rm = T), .groups = "keep") %>%
ungroup() %>%
mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")))# Build Regression Equation Labels
lm_all <- lm(sst_anom ~ year,
data = filter(annual_summary, year %in% c(1982:2020))) %>%
coef() %>%
round(3)
lm_15 <- lm(sst_anom ~ year,
data = filter(annual_summary, year %in% c(2006:2020))) %>%
coef() %>%
round(3)
lm_global <- lm(sst_anom ~ year,
data = filter(global_summary, year %in% c(1982:2020))) %>%
coef() %>%
round(3)
# Convert yearly rate to decadal
decade_all <- lm_all['year'] * 10
decade_15 <- lm_15['year'] * 10
decade_global <- lm_global["year"] * 10
# Equation to paste in
eq_all <- paste0("y = ", lm_all['(Intercept)'], " + ", decade_all, " x")
eq_15 <- paste0("y = ", lm_15['(Intercept)'], " + ", decade_15, " x")
eq_global <- paste0("y = ", lm_15['(Intercept)'], " + ", decade_global, " x")#### Annual Trend Plot ####
ggplot(data = annual_summary, aes(yr_as_dtime, sst_anom)) +
# Add daily data
geom_line(data = region_timeseries,
aes(time, sst_anom),
alpha = 0.5, color = "gray") +
# Overlay yearly means
geom_line(alpha = 0.7) +
geom_point(alpha = 0.7) +
# Add regression lines
geom_smooth(data = filter(global_summary, year <= 2020),
method = "lm",
aes(color = "1982-2020 Global Trend"),
formula = y ~ x, se = F,
linetype = 3) +
geom_smooth(data = filter(annual_summary, year <= 2020),
method = "lm",
aes(color = "1982-2020 Regional Trend"),
formula = y ~ x, se = F,
linetype = 2) +
geom_smooth(data = filter(annual_summary, year %in% c(2006:2020)),
method = "lm",
aes(color = "2006-2020 Regional Trend"),
formula = y ~ x, se = F,
linetype = 2) +
# Manually add equations so they show yearly not daily coeff
geom_text(data = data.frame(),
aes(label = eq_all, x = min(region_timeseries$time), y = Inf),
hjust = 0, vjust = 2, color = gmri_cols("gmri blue")) +
geom_text(data = data.frame(),
aes(label = eq_15, x = min(region_timeseries$time), y = Inf),
hjust = 0, vjust = 3.5, color = gmri_cols("orange")) +
geom_text(data = data.frame(),
aes(label = eq_global, x = min(region_timeseries$time), y = Inf),
hjust = 0, vjust = 5, color = gmri_cols("green")) +
# Colors
scale_color_manual(values = c(
"1982-2020 Regional Trend" = as.character(gmri_cols("gmri blue")),
"2006-2020 Regional Trend" = as.character(gmri_cols("orange")),
"1982-2020 Global Trend" = as.character(gmri_cols("green")))) +
# labels + theme
labs(x = "",
y = expression("Sea Surface Temperature Anomaly"~degree~C),
caption = paste0("Regression coefficients reflect decadal change in sea surface temperature.
Anomalies calculated using 1982-2011 reference period.")) +
theme(legend.title = element_blank(),
legend.position = c(0.85, 0.1),
legend.background = element_rect(fill = "transparent"),
legend.key = element_rect(fill = "transparent", color = "transparent"),
panel.grid = element_blank())
Quarterly
# Repeat for the seasons (equal quarters here*)
quarter_summary <- region_timeseries %>%
mutate(year = year(time),
season = factor(quarter(time, fiscal_start = 1)),
season = fct_recode(season,
c("Jan 1 - March 31" = "1"),
c("Apr 1 - Jun 30" = "2"),
c("Jul 1 - Sep 30" = "3"),
c("Oct 1 - Dec 31" = "4"))) %>%
group_by(year, season) %>%
summarise(sst = mean(sst, na.rm = T),
sst_anom = mean(sst_anom, na.rm = T),
.groups = "keep")
# Plot
quarter_summary %>%
ggplot(aes(year, sst_anom)) +
geom_line(group = 1) +
geom_point() +
geom_smooth(method = "lm",
aes(color = "Regional Trend"),
formula = y ~ x, se = F, linetype = 2) +
stat_poly_eq(formula = y ~ x,
color = gmri_cols("gmri blue"),
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = T) +
scale_color_manual(values = c("Regional Trend" = as.character(gmri_cols("orange")))) +
labs(x = "",
y = expression("Sea Surface Temperature Anomaly"~degree~C),
caption = "Regression coefficients reflect annual change in sea surface temperature.") +
theme(legend.title = element_blank(),
legend.position = c(0.875, 0.05),
legend.background = element_rect(fill = "transparent"),
legend.key = element_rect(fill = "transparent", color = "transparent")) +
facet_wrap(~season)
Marine Heatwaves
For the figures below heatwave events were determined using the methods of Hobday et al. 2016 and implemented using the R package {heatwaveR}
The {heatwaveR} package provides a relatively quick way of working with tabular data to calculate a seasonal climate mean, as well as track heatwave and coldspell events that pass a desired threshold.
These functions were used to track heatwave and cold spell events, counting their frequencies, and tracking their durations. Heatwave events follow the definition of Hobday et al. 2016:
A marine heatwave is defined a when seawater temperatures exceed a seasonally-varying threshold (usually the 90th percentile) for at least 5 consecutive days. Successive heatwaves with gaps of 2 days or less are considered part of the same event.
The heatwave threshold used below was 90%. The heatwave history for Apershing Gulf Of Maine is displayed below:
# Use function to process heatwave data for plotting
region_heatwaves <- pull_heatwave_events(region_timeseries, threshold = 90) %>%
distinct(time, .keep_all = T)Interactive SST Timeline
For anything we wish to host on the web there is an option to display tables and graphs that are interactive. The {plotly} package is one-such tool for producing plots that allow users to pan, zoom, and highlight discrete observations.
# Grab data from the most recent year through present day to plot
last_year <- Sys.Date() - 365
last_year <- last_year - yday(last_year) + 1
last_yr_heatwaves <- region_heatwaves %>%
filter(time >= last_year)
# Get number of heatwave events and total heatwave days for last year
# How many heatwave events:
last_full_yr <- last_yr_heatwaves %>% filter(year(time) == year(last_year))
num_events <- max(last_full_yr$mhw_event_no, na.rm = T) - min(last_full_yr$mhw_event_no, na.rm = T)
# How many heatwave days
num_hw_days <- sum(last_full_yr$mhw_event, na.rm = T)# Plot the interactive timeseries
last_yr_heatwaves %>%
filter(time >= last_year) %>%
plotly_mhw_plots()Static SST Timeline
# Set colors by name
color_vals <- c(
"Sea Surface Temperature" = "royalblue",
"Heatwave Event" = "darkred",
"Cold Spell Event" = "lightblue",
"MHW Threshold" = "coral3",
#"MHW Threshold" = "gray30",
"MCS Threshold" = "skyblue",
#"MCS Threshold" = "gray30",
"Daily Climatology" = "gray30")
# Set the label with degree symbol
ylab <- expression("Sea Surface Temperature"~degree~C)
# Plot the last 365 days
ggplot(last_yr_heatwaves, aes(x = time)) +
geom_segment(aes(x = time, xend = time, y = seas, yend = sst),
color = "royalblue", alpha = 0.25) +
geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe),
color = "darkred", alpha = 0.25) +
geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
geom_line(aes(y = hwe, color = "Heatwave Event")) +
geom_line(aes(y = cse, color = "Cold Spell Event")) +
geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 3, size = .5) +
geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
scale_color_manual(values = color_vals) +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
theme(legend.title = element_blank(),
legend.position = "bottom") +
labs(x = "",
y = ylab,
caption = paste0("X-axis start date: ", last_year,
"\nNumber of Heatwave events in ",
year(last_year),
": ",
num_events,
"\nNumber of Heatwave Days in ",
year(last_year),
": ",
num_hw_days,
"\nClimate reference period : 1982-2011"))
Heatwave Events
# Prep the legend title
guide_lab <- expression("Sea Surface Temperature Anomaly"~degree~C)
# Set new axis dimensions, y = year, x = day within year
# use a flate_date so that they don't stair step
base_date <- as.Date("2000-01-01")
grid_data <- region_heatwaves %>%
mutate(year = year(time),
yday = yday(time),
flat_date = as.Date(yday-1, origin = base_date))
# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(grid_data$sst_anom)) *c(-1,1)
# Assemble heatmap plot
heatwave_heatmap <- ggplot(grid_data, aes(x = flat_date, y = year)) +
# background box fill for missing dates
geom_rect(xmin = base_date, xmax = base_date + 365,
ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5,
fill = "gray75", color = "transparent") +
# tile for sst colors
geom_tile(aes(fill = sst_anom)) +
# points for heatwave events
geom_point(data = filter(grid_data, mhw_event == TRUE),
aes(x = flat_date, y = year), size = .25) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
scale_y_continuous(limits = c(1980.5, 2021.5), expand = c(0,0)) +
labs(x = "",
y = "",
"\nClimate reference period : 1982-2011") +
#scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
scale_fill_distiller(palette = "RdBu", na.value = "transparent",
limit = limit) +
#5 inches is default rmarkdown height for barheight
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "right",
title.hjust = 0.5,
barheight = unit(4.8, "inches"),
frame.colour = "black",
ticks.colour = "black")) +
theme_classic() +
theme(legend.title = element_text(angle = 90))
#### summary side plots:
# number of heatwaves
# remove NA as a distinct heatwave number
n_waves <- grid_data %>%
group_by(year(time), mhw_event_no) %>%
summarise(total_days = sum(mhw_event, na.rm = T), .groups = "keep") %>%
ungroup() %>%
drop_na() %>%
group_by(`year(time)`) %>%
summarise(num_waves = n_distinct(mhw_event_no), .groups = "keep")
# average heatwave duration
wave_len <- grid_data %>%
group_by(year(time), mhw_event_no) %>%
summarise(total_days = sum(mhw_event, na.rm = T), .groups = "keep") %>%
ungroup() %>%
drop_na() %>%
group_by(`year(time)`) %>%
summarise(avg_length = mean(total_days), .groups = "keep")
wave_summary <- left_join(n_waves, wave_len)
# Assemble pieces
heatwave_heatmap
Sea Surface Temperature Anomaly Maps
The 2020 global sea surface temperature anomalies have been loaded and displayed below to visualize how different areas of the ocean experience swings in temperature.
# Access information to netcdf on box
nc_year <- "2020"
anom_path <- str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", nc_year, ".nc")
# Load 2020 as stack
anoms_2020 <- stack(anom_path)Mean Global
The following code will subset the anomalies for the last full month and plot the average sea surface temperature anomalies for that month:
# Get the mean temperature anomalies for the last full month
which_month <- str_sub(Sys.Date(), 6, 7)
which_month <- str_pad(as.numeric(which_month) - 1, width = 2, side = "left", pad = "0")
# Pull last month from the anomaly stack
last_month_dates <- which(str_sub(names(anoms_2020), 7, 8) == which_month)
last_month_avg <- mean(anoms_2020[[last_month_dates]])
# Convert wgs84 raster to stars array
last_month_st <- st_as_stars(rotate(last_month_avg))
# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(last_month_st[[1]]), na.rm = T) * c(-1,1)
# Plot global map
ggplot() +
geom_stars(data = last_month_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_distiller(palette = "RdBu", na.value = "transparent", limit = limit) +
map_theme +
coord_sf(expand = FALSE) +
guides("fill" = guide_colorbar(title = "Mean Sea Surface Temperature Anomaly",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black"))
Regional Map
Using the July average sst anomalies we can then clip the data to just the Apershing Gulf Of Maine for an aerial view of the region. For this execution the bounding box of lat-lon coordinates from the regional extent will be used.
# Clip Raster - Convert to stars
shape_extent <- c(crop_x, crop_y)
region_ras <- crop(rotate(last_month_avg), extent(shape_extent))
region_st <- st_as_stars(region_ras)
# Get crop bounds for coord_sf
crop_x <- st_bbox(region_st)[c(1,3)]
crop_y <- st_bbox(region_st)[c(2,4)]
# Zoom out some for cpr extent to be the same as Andy's GOM
if(tolower(params$region) == "cpr gulf of maine"){
crop_x <- c(-71, -65.5)
crop_y <- c(40.5, 45)}
# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(region_st[[1]]), na.rm = T) * c(-1,1)
# Plot Region - WGS84 Projection
ggplot() +
geom_stars(data = region_st) +
geom_sf(data = new_england, fill = "gray90") +
geom_sf(data = canada, fill = "gray90") +
geom_sf(data = greenland, fill = "gray90") +
scale_fill_distiller(palette = "RdBu", na.value = "transparent", limit = limit) +
map_theme +
coord_sf(xlim = crop_x, ylim = crop_y, expand = T) +
guides("fill" = guide_colorbar(
title = "Average Sea Surface Temperature Anomaly",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black"))
Heatwave Progression
Currently in development. Idea is to animate heatwave through the most recent heatwave event.
# Pull the dates of the most recent heatwave
last_event <- max(region_heatwaves$mhw_event_no, na.rm = T)
last_event_dates <- region_heatwaves %>%
filter(mhw_event_no == last_event) %>%
pull(time)
# Buffer the dates
event_start <- (min(last_event_dates) - 7)
event_stop <- max(last_event_dates)
date_seq <- seq.Date(from = event_start,
to = event_stop,
by = 1)
# Load the heatwave dates
data_window <- data.frame(time = c(min(date_seq) , max(date_seq) ),
lon = crop_x,
lat = crop_y)
# Pull data
hw_stack <- oisst_window_load(oisst_path = oisst_path,
data_window = data_window,
anomalies = T)
#drop any empty years that bug in
hw_stack <- hw_stack[map(hw_stack, class) != "character"]##### Format the layers and loop through the maps ####
# Grab only current year, format dates
this_yr <- stack(hw_stack)
day_count <- length(names(this_yr))
day_labs <- str_replace_all(names(this_yr), "[.]","-")
day_labs <- str_replace_all(day_labs, "X", "")
day_count <- c(1:day_count) %>% setNames(day_labs)
# TESTING: Progress through daily timeline to indicate date
hw_timeline <- cellStats(this_yr, mean)
hw_timeline <- data.frame(date = as.Date(day_labs),
sst_anom = hw_timeline)
# Set palette limits to center it on 0 with scale_fill_distiller
limit <- c(max(values(this_yr), na.rm = T) * -1,
max(values(this_yr), na.rm = T) )
# Plot Heatwave 1 day at a time
day_plots <- imap(day_count, function(date_index, date_label) {
# grab dates
heatwaves_st <- st_as_stars(this_yr[[date_index]])
#### 1. Map the Anomalies in Space
day_plot <- ggplot() +
geom_stars(data = heatwaves_st) +
geom_sf(data = new_england, fill = "gray90") +
geom_sf(data = canada, fill = "gray90") +
geom_sf(data = greenland, fill = "gray90") +
scale_fill_distiller(palette = "RdBu",
na.value = "transparent",
limit = limit) +
map_theme +
coord_sf(xlim = crop_x, ylim = crop_y, expand = T) +
guides("fill" = guide_colorbar(
title = expression("Sea Surface Temperature Anomaly"~degree~C),
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black"))
#### 2. Plot the day and the overall anomaly to track dates
date_timeline <- ggplot() +
geom_line(data = hw_timeline, aes(date, sst_anom),
color = gmri_cols("gmri blue"),
show.legend = F) +
geom_point(data = filter(hw_timeline, date == as.Date(date_label)),
aes(date, sst_anom),
color = gmri_cols("orange"),
size = 3) +
labs(x = "", y = "SST Anomaly")
#### 3. Assemble plot
plot_agg <- (date_timeline / day_plot) + plot_layout(heights = c(1, 3))
return(plot_agg )
})
walk(day_plots, print)
Front Progression
Same idea as above but looking at the fronts rather than absolute values
# remotes::install_github("galuardi/boaR",
# force = T,
# build = T,
# dependencies = "ask",
# upgrade = "ask")
library(boaR)
# front test
front_test <- this_yr[[1]]
front_matrix <- as.matrix(front_test)
dim(front_matrix)
# getting coordinates
test_coords <- raster::coordinates(front_test)
xcoords <- sort(unique(test_coords[,"x"]))
ycoords <- sort(unique(test_coords[,"y"]))
rownames(front_matrix) <- ycoords
colnames(front_matrix) <- xcoords
image(front_matrix)
# Getting Belkin O'Reilly Fronts
boaR::boa(lon = test_coords[,"x"],
lat = test_coords[,"y"],
ingrid = front_matrix,
nodata = NA,
direction = FALSE)
##########################
# Andrew's code for Fronts
## Land and Biophysical regions
shape_path <- shared.path(os.use = "unix", group = "RES_Data", folder = "Shapefiles")
gom_poly<- st_read(paste0(shape_path, "GulfofMainePhysioRegions/PhysioRegions_WGS84.shp"))
focal_area <- st_bbox(gom_poly)
land <- st_read(paste0(shape_path, "ne_50m_land/ne_50m_land.shp"))
# Collecting variables of interest ----------------------------------------
focal_dates<- as.Date(c("2019-04-01", "2019-09-30"))
## Depth - https://www.ngdc.noaa.gov/mgg/global/global.html
depth_info <- rerddap::info("etopo180")
depth <- rxtracto_3D(dataInfo = depth_info,
parameter = "altitude",
xcoord = focal_area[c(1,3)],
ycoord = focal_area[c(2,4)])
depth_rast<- extracto_to_rast(extracto_obj = depth, variable = "depth")
## SST - Daily MUR SST, Final product (MUR-JPL-L4-GLOB-v4.1): https://podaac.jpl.nasa.gov/dataset/MUR-JPL-L4-GLOB-v4.1?ids=&values=&search=MUR (0.01 deg grid, 1 km)
sst_info <- rerddap::info("nasa_jpl_28d8_bd66_24b1")
sst <- rxtracto_3D(dataInfo = sst_info,
parameter = "analysed_sst",
xcoord = as.numeric(focal_area[c(1,3)]),
ycoord = as.numeric(focal_area[c(2,4)]),
tcoord = focal_dates, tName = "time")
sst_rast <- extracto_to_rast(extracto_obj = sst,
variable = "analysed_sst")Warming Rates
Global warming rates are calculated by getting the average temperature for each year across all cells in the data, then using a linear regression for each to determine the annual warming rate in degrees Celsius.
The warming rates are then ranked from high to low to identify areas with the most rapid ocean warming. The following maps display the warming rates and the percentile ranks for the top 20% of warming rates, or rates above the 80th percentile of the data.
# Access information to warming rate netcdf files on box
rates_path <- paste0(oisst_path, "warming_rates/annual_warming_rates")
# Percentile cutoff for displaying warming rates on map
percentile_cutoff <- 801982-2011
#### Data Prep ####
# 1982-2011
rates_stack_all <- stack(str_c(rates_path, "1982to2011.nc"),
varname = "annual_warming_rate")
ranks_stack_all <- stack(str_c(rates_path, "1982to2011.nc"),
varname = "rate_percentile")
# Reclassify to discrete scale
rates_prepped <- reclassify_to_discrete(rates_stack = rates_stack_all,
ranks_stack = ranks_stack_all,
percentile_cutoff = percentile_cutoff)
# pull out components
rates_st <- rates_prepped$rates
rates_raw_st <- rates_prepped$rates_raw
ranks_st <- rates_prepped$ranks
# Get the average warming rate and rank across the area:
rates_crop <- crop(rotate(rates_stack_all), region_extent)
region_rate <- round(cellStats(rates_crop, mean), 3)
ranks_crop <- crop(rotate(ranks_stack_all), region_extent)
region_rank <- round(cellStats(ranks_crop, mean) , 3)The averaged sea surface warming rate percentile for the Apershing Gulf Of Maine over the 1982-2011 period was: 0.624. With a regional warming rate of 0.017 degrees Celsius per year.
Global Rates
#### Plot Setup ####
# guide label
guide_lab <- expression("Annual Warming Rate - "~degree~C)
# ggplot using stars
ggplot() +
geom_stars(data = rates_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
"\n Data from Years: 1982-2020"))
Global Ranks
# ggplot using stars
ggplot() +
geom_stars(data = ranks_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = "Percentile of Annual Warming Rate",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
"\n Data from Years: 1982-2020"))
2006-2020
#### Data Prep ####
# 2006-2020
rates_stack_all <- stack(str_c(rates_path, "2006to2020.nc"),
varname = "annual_warming_rate")
ranks_stack_all <- stack(str_c(rates_path, "2006to2020.nc"),
varname = "rate_percentile")
rates_prepped <- reclassify_to_discrete(rates_stack = rates_stack_all,
ranks_stack = ranks_stack_all)
rates_st <- rates_prepped$rates
ranks_st <- rates_prepped$ranks
# Get the average warming rate and rank across the area:
rates_crop <- crop(rotate(rates_stack_all), region_extent)
region_rate <- round(cellStats(rates_crop, mean), 3)
ranks_crop <- crop(rotate(ranks_stack_all), region_extent)
region_rank <- round(cellStats(ranks_crop, mean) , 3)The averaged sea surface warming rate percentile for the Apershing Gulf Of Maine over the 2006-2020 period was: 0.916. With a regional warming rate of 0.077 degrees Celsius per year.
Global Rates
#### Plot Setup ####
# guide label
guide_lab <- expression("Annual Warming Rate - "~degree~C)
# ggplot using stars
ggplot() +
geom_stars(data = rates_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
"\n Data from Years: 2006-2020"))
Global Ranks
# ggplot using stars
ggplot() +
geom_stars(data = ranks_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = "Percentile of Annual Warming Rate",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
"\n Data from Years: 2006-2020"))
1982-2020
#### Data Prep ####
# 1982-2020
rates_stack_all <- stack(str_c(rates_path, "1982to2020.nc"),
varname = "annual_warming_rate")
ranks_stack_all <- stack(str_c(rates_path, "1982to2020.nc"),
varname = "rate_percentile")
rates_prepped <- reclassify_to_discrete(rates_stack = rates_stack_all,
ranks_stack = ranks_stack_all,
percentile_cutoff = percentile_cutoff)
rates_st <- rates_prepped$rates
rates_raw_st <- rates_prepped$rates_raw
ranks_st <- rates_prepped$ranks
# Get the average warming rate and rank across the area:
rates_crop <- crop(rotate(rates_stack_all), region_extent)
region_rate <- round(cellStats(rates_crop, mean), 3)
ranks_crop <- crop(rotate(ranks_stack_all), region_extent)
region_rank <- round(cellStats(ranks_crop, mean) , 3)The averaged sea surface warming rate percentile for the Apershing Gulf Of Maine over the 1982-2020 period was: 0.944. With a regional warming rate of 0.041 degrees Celsius per year.
Global Rates
#### Plot Setup ####
# guide label
guide_lab <- expression("Annual Warming Rate - "~degree~C)
# ggplot using stars
ggplot() +
geom_stars(data = rates_st) +
#geom_stars(data = rates_raw_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
"\n Data from Years: 1982-2020"))
Global Ranks
# ggplot using stars
ggplot() +
geom_stars(data = ranks_st) +
geom_sf(data = world, fill = "gray90") +
scale_fill_viridis_c(na.value = "black", option = "plasma") +
guides("fill" = guide_colorbar(title = "Percentile of Annual Warming Rate",
title.position = "top",
title.hjust = 0.5,
barwidth = unit(4, "in"),
frame.colour = "black",
ticks.colour = "black")) +
map_theme +
coord_sf(expand = FALSE) +
labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
"\n Data from Years: 1982-2020"))
A work by Adam A. Kemberling
Akemberling@gmri.org